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摘要 经 典 的 初 轨 确 定 方 法 包括 Laplace 方 法 和 Gauss 方 法 以 及 它们 的 各 种 变化 形式 . 除 这 些 经 典 方法 之 外 ， 
| 于 当今 光学 观测 数据 的 特点 ， 学 者 们 也 陆续 提出 了 一 些 其 他 的 初 轨 确 定 方法 , 包括 双 r (目标 距离 观测 者 的 距 
离 ) 方 法 和 可 行 域 方 法 . 双 7 方 法 的 一 种 实现 方式 是 通过 猜测 某 两 个 时 刻 (通常 是 定 轨 弧 段 的 首 、 末 时 刻 ) 目 标 离 
观测 者 的 距离 , 结合 观测 者 在 空间 中 的 位 置 矢量 , 即 可 求解 相应 的 Lambert 弧 段 作为 目标 轨道 的 初始 猜测 . 进 
步 , 以 其 他 观测 时 刻 的 RMS (Root Mean Square) 为 优化 变量 可 以 改进 初始 猜测 从 而 确定 初 轨 . 可 行 域 方法 则 是 
针对 一 组 初始 观测 参数 (包括 赤 经 、 赤 纬 及 其 变 率 ), 根据 一 些 初始 假设 将 目标 ( 离 观测 者 的 ) 距 离 及 其 变 率 约束 
在 可 行 域内 , 并 通过 三 角 划 分 逐步 逼近 的 方式 寻找 到 使 观测 RMS 最 小 的 猜测 解 . 针对 一 系列 模拟 观测 数据 以 及 
实测 数据 , 将 智能 优化 算法 (粒子 群 算法 ) 应 用 于 这 两 种 初 轨 方 法 , 并 将 结果 与 改进 的 Laplace 算 法 的 结果 进行 比 
BE. 由 于 双 7 方 法 不 仅 可 以 用 于 短 弧 定 轨 还 可 用 于 长 弧 关联 , 所 以 进一步 给 出 了 针对 长 弧 段 数据 的 关联 结果 . 
关键 词 天 体力 学 , 小 天 体 , 方法 : 数值 
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1 5| 言 工作 中 , 我 们 只 考虑 基于 近 地 天 体 光学 观测 资料 的 
通常 来 说 , 近 地 天 体 的 轨道 确定 包括 初 轨 确 定 IEA. N 

= 和 精密 定 轨 两 个 步骤 , 初 轨 确 定 是 利用 较 少 的 光学 从 1780 年 Laplace 方 法 被 提出 至 今 ， 经 典 的 

= 或 雷达 观测 资料 , 为 精密 定 轨 提供 轨道 初 值 . 初 轨 初 轨 计算 方法 主要 有 Laplace、 Gauss. Escobal 

= 确定 的 结果 越 好 , 精密 定 轨 的 过 程 越 容易 收敛 . — IK, "是 目标 距离 观测 者 的 距离 ) 方 法 以 及 它们 


般 来 说 , 光学 观测 资料 包括 目标 的 赤 经 和 赤 纬 ,而 的 衍生 方法 . Laplace 方 法 和 Gauss 方 法 都 是 对 3 组 
随 着 观测 技术 的 发 展 , 出 现 了 光子 计数 的 探测 装置 ， 光学 观测 资料 进行 公式 推导 ， ih Gace 


通过 这 种 装置 得 到 的 原始 数据 , 除 角 度 信 息 外 还 可 ”一 个 八 阶 方程 的 求解 . PRT, Charlier!” 
获得 角度 变化 率 信息 , 在 很 多 情况 下 , 即使 是 短 弧 “指出 Laplace 定 初 轨 的 方法 给 出 的 解 有 时 并 不 唯 
段 也 能 得 到 角度 的 二 阶 变 率 信息 凹 . 部 分 近 地 天 体 一 , Gronchi 指 出 此 问题 对 Gauss 定 初 轨 方 法 同样 存 


还 有 雷达 观测 资料 (时 延 和 多 普 勒 频 移 ), 但 在 本 文 EB, 一 种 改进 Laplace 方 法 向 避 开 了 对 八 阶 方程 的 
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求解 ， 而 是 采用 数值 迭代 的 方法 求解 “线性 ”表达 
ES 式 的 非 线性 方程 组 . 1965 年 Escobal 回 提出 双 r 方 
YE, 之 后 Goodingl 四 和 Briggs 等 门 分 别 对 双 r 方 法 进 
行 了 改进 . 双 r 方 法 的 本 质 都 是 假设 已 知 两 组 观 
测 资料 的 径 向 信息 (目标 与 观测 者 的 距离 ) 即 可 
根据 其 余 的 观测 资料 对 这 两 组 径 向 信息 进行 修 
IE; "提出 的 方法 借助 了 轨道 积分 方程 ， 
而 Goodingls| 借 助 了 求解 Lambert 弧 段 的 方法 . 除 
此 之 外 , 还 有 多 种 不 同情 况 下 不 同方 式 的 改进 方法 ， 
壁 如 Zeinalov 四 在 只 考虑 圆 轨道 的 情况 下 提出 了 没 
9 时间 信息 情况 下 的 解决 方案 . 

近 现 代 由 于 观测 技术 提高 , 出 现 了 大 量 没有 归 


Briggs 等 
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过 求解 兰 伯 特 问题 来 确定 初始 时 刻 的 轨道 参数 L171. 
Hinagawa 等 人 采用 遗传 算法 利用 地 基 光 学 系统 对 
短 弧 段 内 的 地 球 同步 轨道 (Geostationary Orbit, 
GEO) 目 标 进 行 轨道 确定 , 将 轨道 根 数 分 为 两 组 作 
为 优化 变量 分 别 进行 求解 (站 , ERE HE MR T 
遗传 算法 优化 求解 了 空间 目标 的 极 短 弧 定 轨 问 题 ， 
这 些 工作 同样 选择 了 轨道 根 数 而 非 视 向 距离 作为 
MBB 271 此 外 , 李 和 独自 等 人 还 针对 近 地 小 行 
E, 采用 进化 算法 进行 了 初 轨 计 算 , 同样 也 是 选择 
初始 轨道 根 数 作为 优化 参数 , 并 提出 对 于 偏心 率 大 
的 天 体 需 要 一 定 程度 分 区 域 寻找 最 优 解 B39. 

本 文 将 考察 近 地 小 天 体 的 初 轨 确 定 ， 针 对 
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类 的 极 短 弧 段 , 此 时 将 经 典 的 定 初 轨 方法 直接 应 用 
于 单个 极 短 弧 段 通 常 难以 奏效 , 不 能 获得 有 效 的 轨 
道 信 息 . 基于 当今 光学 观测 资料 的 特点 (有 较为 准 
确 的 角度 变 率 信息 )，Milani 等 鲁 针 对 小 天 体 提出 了 
可 行 域 方法 , 此 后 Tommei 等 人 又 将 此 方法 应 用 到 
近 地 空 间 目标 的 初 轨 确定 问题 , 并 推广 到 雷达 类 型 
观测 资料 09. 以 光学 资料 为 例 , 这 种 方法 采用 的 基 
本 观测 集 (TSA 一 Too Short Arc) 除 了 两 个 角度 之 
外 还 包含 它们 的 变 率 . 对 一 个 TSA, 根据 一 定 约束 
段 设 可 以 确定 可 行 域 . 通过 对 两 个 TSA 的 可 行 域 进 
行 关 联 可 以 进一步 确定 一 条 初始 轨道 器. 在 不 同 应 
情况 下 , 学 者 们 给 出 了 不 同 的 限制 , 比如 Marusk- 
in 等 [限制 了 近地点 和 远地点 ，DeMars 等 (站 又 根 
据 具体 情况 给 出 了 半 长 径 和 偏心 率 的 限制 . 值得 
提 的 是 , 另 一 种 基于 两 个 TSA 的 初 轨 确定 方案 是 直 
接 借助 能 量 积 分 和 角 动 量 积 分 求解 上 LI . 

智能 计算 也 有 人 称 之 为 “ 软 计算 ”, 是 人 们 受 自 
然 ( 生 物 界 ) 规 律 的 启迪 , 根据 其 原理 , 模仿 求解 问 
题 的 算法 . 智能 优化 算法 有 效 解决 了 在 多 变量 多 
约束 情况 下 的 寻 优 问题 ,常用 的 有 退火 算法 、 遗 
传 算 法 、 粒 子 群 算法 等 . 有 关 学 者 已 将 智能 优化 
算法 融入 了 初 轨 计算 中 . 王 雪 莹 等 人 将 粒子 群 优 
化 算法 和 可 行 区 域 相 结合 ,不 再 将 可 行 区 域 进行 
三 角 化 , 而 是 直接 利用 粒子 群 算法 对 天 基 短 弧 观 
测 数据 进行 全 局 寻 优 , 提高 了 计算 效率 l16]. Ansa- 
lone 等 人 将 遗传 算法 和 双 r 方 法 相 结合 , 对 极 短 弧 段 
天 基 观 测 数据 下 的 空间 目标 轨道 参数 进行 估计 , 将 
观测 弧 段 内 首 末 时 刻 的 视 向 距离 作为 优化 参数 , 通 


44-2 


偏心 率 较 大 ( 约 为 0.6) 的 近 地 小 天 体 Toutatis ( 编 
号 4179), 采用 两 种 方法 一 可 行 域 方法 以 及 双 r 方 法 ， 
利用 粒子 群 算法 针对 优化 变量 进行 全 局 寻 优 , 同时 
给 出 了 协 方差 矩阵 的 估算 方法 . 考虑 到 在 不 同 场景 
下 采用 的 弧 段 情况 不 同 , 本 文 还 分 别 计算 不 同时 ] 
和 观测 数目 的 弧 段 , 比较 了 改进 Laplace 方 法 、 可 
域 方 法 以 及 双 r 方 法 的 优 劣 性 ,分析 了 不 同 弧 段 
长 和 观测 数目 对 不 同方 法 精度 的 影响 , 试验 结果 : 
示 弧 段 时 长 对 于 改进 Laplace 方 法 和 双 7 方 法 影 
较 大 , 而 观测 数目 对 可 行 域 方法 的 影响 较 大 . 以 出 
为 基础 , 本 文 最 后 给 出 针对 不 同 长 度 弧 段 进行 初 旬 
计算 所 应 采用 方法 的 建议 . 

本 文 在 第 2 节 简 述 了 从 经 典 到 现代 的 初 轨 计 
算 方 法 一 高 斯 方法 、 改 进 Laplace 方 法 、 双 7 方法 
以 及 可 行 域 方 法 .第 3 节 描 述 了 与 智能 优化 算法 
结合 的 双 7 方 法 和 可 行 域 方 法 ， 并 给 出 包括 改进 
Laplace 方 法 在 内 的 3 种 方法 的 协 方差 矩阵 的 计算 . 
在 第 4 节 , 展示 了 本 文 所 用 的 数据 的 处 理 方式 以 及 
采用 3 种 方法 针对 不 同 弧 段 进 行 初 轨 计 算 的 结果 . 
最 后 对 试验 结果 进行 了 讨论 和 总 结 . 


2 ” 初 轨 方 法 简介 
经 典 初 轨 确 定 方法 
2.1.1 ”Gauss 方 法 

图 1 是 中 心 天 体 太 阳 S、 地 球 五 以 及 小 天 体 4 的 
相对 几何 构 型 ， 其 中 , 7 表示 小 天 体 相 对 于 日 心 的 
MEARE, p 表 示 小 天 体 相对 于 测 站 的 位 置 矢 量 ; 
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R= (X,Y, 2)T 表 示 测 站 相对 于 日 心 的 位 置 矢量 ; 2q° (Ceosez + y)r3 — g? , 


GERRENA ANRA; BARRERA yeh py ERATARA, CARTE GLE 
的 夹 角 . 在 国际 天 球 参考 系 (International Celestial 中 的 常数 , g = Ro oHa SRZ Hkt: 并 且 有 
Reference System, ICRS) 下 的 测量 几何 满足 如 下 关系 : | | 


关系 : 4 
By ar = ay, 
s C C 3 
F=p+R. (1) ò J (3) 
C= 一 1 一 3° 
q 72 


(3) 式 又 称 为 动力 学 方程 (Dynamical Equation). 这 
里 假设 y = 1, 有 


: L P(r) = (r2 — a) Para). 


Pi (r2) = Cr$ (r2 +q) + (73 +qr2 +°) 


= r < lg” 一 (2Ccoses + 1)q?r9] . 

oe $ 此 时 

© § PO) == <0). lim P(r2) = +00. 

rT 

= 图 1 太阳 、 地 球 以 及 小 天 体 的 相对 几何 构 型 fea on Ee 可 以 知道 该 方程 一 定 
Fig.1 The relative geometric configuration of the Sun, 有 3 个 正 根 ， 由 于 物 理 性 M, 72 = = 4 的 根 应 被 去 掉 . 
= Earth and Asteroid 另外 两 个 根 根据 该 表达 式 ， 由 零 圈 (zero circle) fll 限 


制 曲 线 (limiting curve) 将 平面 分 为 几 个 部 分 . 图 2 展 

考虑 3 组 观测 资料 {taa 51); (oa 92), (03,63)} 示 了 零 圈 和 限制 曲线 将 整个 平面 划分 不 同 解 的 区 

对 应 观测 里 刻 (二 123). 其 中 Qs、5 分 别 表示 m, 零 圈 以 外 的 区 域 表 示 有 两 个 根 ， 零 圈 和 限制 

第 ;组 观测 数据 的 亦 经 和 亦 纬 , 并 有 对 应 测 站 在 太 。 线 之 间 的 区 域 表 示 一 个 根 , 限制 曲线 以 内 表示 两 个 

阳 系 日 心 天 球 坐标 系 中 的 位 置 矢量 尼 (i = 12,3), 。 根 . 该 部 分 Gronchi 等 四 给 出 了 具体 推导 , 这 里 仅仅 

以 Ri (05 1, 2,3) AER, 多 二 (一 十 23) 表示 小 利用 该 图 说 明 求解 初 负 问题 时 会 出 现 多 解 的 情况 
天 体 在 日 心 天 球 坐 标 系 中 的 位 置 矢量 .Gronchi 竺 
a ian 2.1.2 ”改进 Laplace 方 法 

人 给 出 了 高 斯 方法 的 具体 推导 [9, 此 处 不 再 效 述 , A E E TTEA 

P 不 考虑 摄 动 的 情况 下 , 小 天 体 的 运动 方程 满足 


s__ Bp 
— [4} + 2.42 ( 3 - Ra)/p2 + Ralr3— | oo wt (4) 
2Bo|A2 + (p2 R2)/p2)r3 一 = B? =0, (2) to : rò = r(to), rò = 7(to), 

eae l EHu = GM 表 示 关 于 中 心 天 体 的 引力 常数 ，G 为 
其 中 , 4。 Bs 是 由 已 知 条 件 组 合 得 到 的 系数 ; rK 万 有 引力 常数 ，M 为 中 心 天 体质 量 , R = (zo, Yo, 
示 第 2 个 观测 时 刻 小 天 体位 置 和 撩 量 避 的 模 长 ; pot wn) Hy EL RS CATIA 0 RL, RR 
示 第 2 个 观测 时 刻 由 观测 站 到 观测 目标 的 相对 位 置 。” 矢量 7# 的 两 阶 导数 , 即 该 目标 天 体 的 加 速度 , F(t) 表 
RE. 利用 牛顿 迭代 的 方法 即 可 得 到 该 方程 的 解 。 ” 示 在 t 时 刻 目标 天 体 的 速度 . 在 短 弧 的 前 提 下 , 该 函 


由 于 Laplace 方 法 最 终 也 会 得 到 形 如 (2) 式 所 示 的 八 数 满足 初始 条 件 的 解 可 以 展开 成 时 间 间 隔 At = t 
阶 方程 , (2) 式 通常 又 写成 下 面 更 普遍 的 形式 一 如 的 震级 解 
1 
P(r2) = C?r§ — (C? + 2Cycosez + 7° )r$+ F(t) = rò + TOM At + xr AP 十 … :十 
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1 
SmE Ag a... 
g? At + 5 


(5) 


Hepes yr (tee EAE AMON, Bll 
Bee a 2) #3 


由 (4) 式 中 运动 方程 给 出 , 即 
k k > + 
AM 二 FEL (to, Fo, Fo) . 
2 
l 
1 solution 
= zero curve 2 solutions 1 solution 
= 0 |. ------------|-----+------ @.---------- ef --------- 
> uae 
limiting curve 
-1| 2 solutions 
-2 
-2 -1 0 1 
x/AU 
2 ” 零 圈 、 限 制 曲线 划分 解 的 个 数 示意 图 


Fig.2 Schematic diagram of the number of solutions for the 


division of zero circle and limiting curve 


同时 有 动 | = ry, 又 可 以 将 其 展开 为 
F(t) = F* (Fo, Fo, At) + G* (Fo, ro, At)ro, (6) 
F PFG HY DA GE E (6) RM5 (5) RA E GE a 
到 , 具体 形式 如 下 : 
1 1 
Fe =1- zur + 本 wopor + 
1 1 5 
(rom T =r hte, 
1 1 
Gr =T- guor + 了 Twopor free, 
(7) 


其 中 wo = 1/r3. po = (ro - 7o)/72. do = 03/3, fo 
表示 而 的 模 长 , 上 述 参 数 下 标 0 都 表示 初始 时 刻 to 时 
的 值 . 
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将 测量 几何 关系 (1) 式 和 (6) 式 结合 并 同时 又 乘 

5 的 单位 矢量 5 可 得 


一 


px [F* (Fo, Fo, At)Fo + G” (Fo, To, At)ro| = px R. 


对 于 3 组 观测 资料 (ai, 6;), 有 几何 关系 
COS Qi COS 6; ri 
Pi = pi | sina; cos 6; =pil m |= pipi 
sin 0; Vi 
由 每 组 观测 资料 对 应 的 测 站 信息 B= (Xa, Yi, Zi) 
可 以 得 到 
F*vi)xo — (F*Ai)zo + (Gšvi)to— 
G% Ài Zo = Vi Xi 3 Zi 3 
(8) 


T 


CIOE aA E AATE 
以 看 到 上 式 中 只 有 两 个 式 子 是 独立 的 , ERARA 
i, 需要 至 少 3 组 观测 资料 , 这 里 记 (8) 式 计算 出 来 的 
ARA E). 利用 (8) 式 不 断 使 用 下 面 公式 进行 
迭代 


) 表 示 第 mm 次 迭代 计算 的 值 


Kp” a 
了 归 一 化 , E 


人 一 VF 


直到 e 满 足 一 定 精度 , 迭代 结束 . 在 实际 工作 中 , 该 
方法 可 对 观测 数据 大 于 等 于 3 组 的 观测 弧 段 有 效 . 


2.2 MrAK 

测 角 资料 确定 初 轨 的 本 质 困难 在 于 缺乏 径 向 
i JS (BES AUK FE), aE i 思想 在 于 知道 
至 少 两 个 观测 点 的 径 向 信息 , 即 可 确定 轨道 . 这 利 
方法 既 可 以 用 于 极 短 弧 (Too Short Arc) 定 轨 , 也 可 
用 于 间隔 时 间 很 长 的 (甚至 多 圈 ) 弧 段 资料 的 初 轨 确 


m—1 A {m—1 
本 


Ne 


z; 


=H 
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定 . 通常 这 种 方式 需要 至 少 3 个 点 的 测 角 资料 , 并 一 已 代表 猜测 的 位 置 , 此 时 定义 一 个 平面 , AEP, HE 
般 假设 第 1 个 和 最 后 一 个 测 角 资料 的 径 向 距离 ， 以 直 于 to 时 刻 目标 位 置 的 真实 方向 , 再 以 该 平面 与 该 
此 推导 中 间 时 刻 观测 点 的 径 向 距离 ， 并 作为 基础 方向 的 交点 为 原点 建立 坐标 系 . 令 忆 在 该 坐标 系 上 
进行 修正 . 总 的 来 说 , 双 r 方 法 有 两 种 形式 , 一 种 与 的 投影 为 1 和 g, 坐标 系 的 建立 需 使 得 六 > OFF Ag = 


Lamberti BHAA, 一 种 与 轨道 运动 积分 相 结合 ， 0. 根据 该 坐标 系 的 建立 可 得 , 当 已 是 目标 的 真实 位 
下 面 分 别 对 两 种 方法 进行 简 述 . 置 时 f = 0. Gooding! 3] Newton-RaphsonZ xt 
2.2.1 与 轨道 积分 结合 的 双 7 方 法 进行 迭代 修正 . KH Aa. yt ei. ps, 

这 里 将 简 述 Briggs 等 四 提出 的 方法 ， 该 方法 求 af af\ | 
解 所 谓 的 Gibbs 问 题 , 又 被 称 为 BSE (Binary-Single (x)-- ðs Dy (7 mn 
Evolution) 方 法 . 首先 假设 有 3 组 观测 资料 , 根据 小 Ay og 9g g) 
Ae AIEEE, (i = 1,2,3) 共 面 的 性 质 以 Ox Oy 
及 轨道 方程 , 可 以 推导 得 到 6 个 轨道 根 数 的 表达 式 ， 其 中 f、g 对 z、y 的 偏 导数 通过 数值 差分 方法 给 出 . 


在 此 不 进行 次 述 , 而 平 近 点 角 MMi; (i = 1,2, 3) 满 足 Gooding9 提 供 的 方法 比 Briggs 等 四 提供 的 方法 有 
两 个 优点 , 首先 Gooding 上 的 方法 收敛 域 更 广 , 采用 


> Aa 求解 Lambert 弧 段 的 方法 不 管 提供 怎样 的 初 值 都 能 
© ” 其 中 m 表 示 平 运动 速度 , w 表 示 近 日 点 辐 角 . 令 w = ”得 到 一 个 合适 的 解 ; 其 次 , 如 果 本 身 的 解 不 确定 , 经 
一 nT, Wt; = —T, +“, 即 可 计算 得 到 t13、tos. © 典 的 方法 与 BSE 方 法 都 无 法 收敛 , 本 文 将 对 这 种 方 
Oo x 法 加 以 简单 改进 以 适用 于 智能 优化 算法 
中 与 观测 的 时 间 可 以 作 差 . 于 是 该 问题 转化 成 求 p1 和 e : 
和 六 使 得 目标 函数 满足 pr 
gh fin fey = 0; k=1,2. (9) i i. 
c 
上 式 可 采用 Newton-Raphson 公 式 进 行 迭 代 求 解 . if \ ` 
Ai(p1, Pp2)、A2(p1,p2) 描 述 了 pi 一 p2 平 面 上 的 两 条 # | 


O WA, 而 这 两 条 曲线 相交 的 点 ， 即 为 该 问题 的 解 . Fa Ps 
a Briggs NEXE P HR T RL A AR d 
== 结束 条 件 , 同时 也 提出 了 在 受到 扰动 的 情况 下 的 修 图 3 Gooding[ 采 用 坐标 系 的 示意 图 . 星 号 表示 PP. 的 位 置 ,虚线 表示 
oad E, XARA. 第 2 个 观测 时 刻 的 方向 . PUMA eae P. ELE TARE. 箭头 
os Statin dee ae 表示 在 新 平面 上 所 建立 的 坐标 系 . 
这 里 将 基于 Goodingls| 提 出 的 方法 进行 简 述 . 


D! 
Ps 


Fig.3 Schematic diagram of the coordinate system used by 


Goodingl6]. The star represents the position of 已 . The 


简单 来 说 ， Gibbs 问题 即 通过 迭代 求解 m 和 po 的 值 dashed line represents the direction at the second 
满足 约束 条 件 ; 而 在 Goodingl6l 提出 的 算法 中 ， 仍 observation time. The shaded area represents the plane 
bp LG > 求解 两 个 占 的 径 Br i k te yt passing through P. and perpendicular to the dashed line. 
pak FTE ACR AE PA | PARI fr 向 距 pit P3 的 (EL The arrow represents the coordinate system established on 
足 两 个 约束 条 件 ， 但 它们 的 定义 不 同 . the new plane. 


在 猜测 了 pi 和 ps 的 值 后 ， 可 以 利用 Lambert 问 
题 求解 该 轨道 , 再 根据 中 间 时 刻 可 得 到 ts 时 刻 的 位 2.3 ”可 行 域 方法 
AP. 图 3 中 虚 曲 线 代表 真 实 轨道 Pi. Py PAR 得 益 于 观测 设备 的 改进 和 较 大 规模 的 巡天 计 
表 目 标 真 实 位 置 , MAAN AIA AIA. Pi 划 , 现 有 的 光学 观测 模式 已 经 发 生 改 变 . 大 多 数 时 
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候 望 远 镜 指向 天 空 的 固定 区 域 短 时 曝光 几 幅 图 片 ， 4. 与 地 球 距 离 大 于 地 球 半 径 


然后 再 指向 不 同 的 天 区 进行 观测 . 近 地 天 体 的 发 
现 通常 采用 对 比 同一 天 区 的 几 幅 照片 的 方式 , 但 
一 般 而 言 , 新 的 大 型 巡天 设备 无 法 对 目标 进行 后 续 其 中 产 ? 瑚 示 太 阳 系 日 心 坐 标 系 下 小 天 体 的 位 
的 跟踪 观测 ,只 能 获得 一 晚 的 数据 P3, 而 这 么 短 的 速度 矢量 ，Es 表示 小 天 体 在 该 坐标 系 下 的 能 

时 长 不 足以 利用 经 典 的 初 轨 确定 方法 进行 初 轨 确 元 、 售 表示 地 心 天 球 坐 标 系 下 小 天 体 的 位 置 、i 
E. 这 么 短 的 观测 弧 段 ,也 被 称 为 极 短 弧 . 对 于 短 ERE, Mre ret am LBW eB, 
时 曝光 的 几 幅 图 片 , 通常 根据 几 幅 照片 的 插值 获得 Erasty 表 示 小 天 体 在 该 坐标 系 下 的 能 量 . Rs 表示 地 
某 个 参考 时 刻下 对 应 的 (a, 6, G6)， 即 运动 目标 在 。 球 半径 , ue = GMe, us = GMs，Ms、Ms 分 别 
站 心 天 球 坐 标 系 下 的 赤 经 和 赤 纬 以 及 他 们 的 变 率 ， 表示 地 球 、 太 阳 质 量 , ae 表示 地 球 公转 轨道 的 半 长 
TSA 对 应 的 数据 又 被 称 为 attributable， Milani 等 加 径 . 故 可 行 区 域 为 

提出 了 此 名 称 , 用 来 表示 此 TSA 是 否 可 与 已 知 弧 段 


0 


Rp = { (F, P)|re > Rg}, 


z 


H 


关联 , 并 以 此 为 基础 分 析 了 一 段 TSA (或 者 说 attri- (Re (Ro) J Re) N Raf Ro. 
putable) 中 所 包含 的 信息 ， 并 提出 了 可 行 域 的 要 最 后 将 可 行 区 域 三 角 化 , 每 个 三 角 区 域 的 顶点 
AON, 此 处 将 进行 简单 介绍 . 针对 于 TSA 中 所 包含 “都 能 给 出 (a, d,à, ò, p, A) AE, 如 果 近 地 小 天 体 在 
的 信息 , 这 些 信息 在 日 心 坐 标 系 下 有 如 下 关系 该 项 点， 则 由 此 信息 可 以 给 出 近 地 小 天 体 的 状态 
量 和 轨道 根 数 . 如 果 有 该 小 天 体 的 另 一 个 TSA 的 数 
r=ptR, 据 , 则 可 进行 关联 以 生成 初始 轨道 的 估计 Ga 


p= pp + (Pod SF Pab) ) 


3 结合 智能 优化 算法 的 改进 


一 Sinacos0 —cosasind 
> Be _ 3.1 Gi ge Cine 
Pa = cos a cos 6 , Pe= | —sinasind 
0 ts 假设 有 NR (NR > 3) 组 观测 资料 (ae ,Bo ),， m 
(11) = 1,2,... ,NR, 假定 已 知 
第 1 个 和 第 NR 个 观测 资料 的 站 心 距离 pl 和 pwr， 即 


若 要 得 到 观测 目标 的 轨道 信息 , 还 需要 得 到 径 
向 信息 , 即 测 距 p 和 测速 p. 可 行 域 方法 则 是 利用 目 
标的 物理 信息 进行 区 域 划分 , 下 面 根据 近 地 小 天 体 


可 得 到 小 天 体 的 日 心 位 置 矢量 ( 见 图 1), 由 此 可 求解 
兰 伯 特 问题 并 计算 中 间 各 观测 时 刻 的 理论 观测 值 


E 7 (ac BS) m = 1,2,.… ,NR, 其 中 上 标 c 表 示 计 算 
性质, 给 5 近 地 小 天 体 不 会 
出 可 行 区 域 , 并 假设 近 地 小 天 体 不 会 被 en se 


j Ym = on 2 ae j ô — 6° a 
1. 小 天 体 属 于 太阳 系 Ym = [cos dm (am 一 amh) (Om — 6m)] 


~ 


l wi 所 有 残 差 向 量 记 为 
Ra= 六 su 一 一 一 十 二 7 O}; 
A {(7 r)| S r 9” < } Y = (Ji, J2, INR)". 
2. 小 天 体 不 会 被 地 球 俘获 记 变 量 系统 为 筷 , 在 本 问题 中 
= T 
Rg = {(r, r) |Etartn = -2 + are? > 01; X = (pi, pwr) - 
| 采用 双 r 算 法 的 初 轨 确 定 问题 即 为 寻找 合适 的 
3. 与 地 球 距离 超过 地 球 希 尔 范围 XIE FRKE RME 
| a aaea 
Ro = {(F, Pre > ap(ug/3)"}; min J(X) = zR Y ri (12) 
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相当 于 求解 下 式 
OJ 
引入 和 矩阵 
B (bı, ba, byr) ; (14) 
其 中 ， 
OYm 
bm aX 
(13) 式 可 以 简 记 为 
OJ TD 
Ta =—aY B=0. (15) 
记 (15) 式 的 解 为 X "引入 
加 » PJ T j 
AX=X-X*, ae =a (B B+Y SG 


个 合 


合理 的 假设 是 ， 当 满足 (15) 式 时 ， 观 测 残 差 
Y 一 0, 所 以 


2 
ax? NR 
考虑 此 近似 , (15) 式 的 迭代 公式 变 为 


B'B. 


JIN / ðJ N” 
0 (Se) (sx) = BB Ms 
(16) 
引入 协 方差 矩阵 
Py, 0 0 
0 P Beate 0 
Py =E[YY"]=]| . us l . l 
0 0 Pyra 
其 中 , 卫 表 示 求 数学 期 望 算 符 . WU EX A 
差 矩阵 为 
Py =(B'B)'B'PyB(B'B)*. (17) 


考虑 到 NR 值 可 能 很 大 , 实际 计算 时 我 们 可 以 


采用 


NR-1 
B'B = > b? bm, BUY 


能 优化 算法 的 小 天 体 初 轨 确 定 


NR-1 
= 》 biyn, B’PyB 


m=2 
NR-1 


= So bn Py. bn 


m=2 


避免 大 矩阵 的 直接 计算 . 在 计算 Py 时 除 需 考虑 观 


测量 的 方差 2 外 还 需 考虑 来 自 第 1 个 和 第 NR 个 观 
测 时 刻 的 角度 误差 贡献 , 即 
Py =P +Py, (18) 
其 中 
po — cos? 09,a2 0 
Ym 0 02 


为 观测 量 的 方差 , cs。 、cav LA AAA 
的 方差 . 而 


Pi. = 
o 0 0 0 
dy, 0 of 0 0 
O(a? of an on | 0 0 oa 0 
0 0 0 只 
3Y m g 
Lô(a?, ôT, QRR, ORR) 
为 第 1 个 时 刻 和 第 NR 个 时 刻 的 观测 量 的 方差 映射 
至 中 间 各 观测 时 刻 的 残 差 . 但 是 (17) 式 给 出 的 仅仅 


日 


是 待 估量 兰 的 方差 , 如 需 给 出 初始 时 刻 的 位 置 和 速 
度 对 应 的 协 方差 矩阵 已 ;, 则 需要 进一步 计算 . 记 


X = (X, 04,67, oR, NR) 
= (P1, PNR: 04,51, ANR ONR) > 
则 
Py 0 0 0 ù 
0 2 0 0 0 
Py=| 0 0 of 0 0 
0 0 mi -0 
0 0 0 0 of, 
5] AHEM 
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i 
COS COSQ1 COSONRCOSQNR —/1 COS ol SİN Q1 
cos ol sin al cosdnRsinanR pı COS 01 cosa, 
sin ôi sin ONR 0 
Oz/Op1 Oz /OpNR 0%/0a4 
Oy /Op1 Oy/Opnr dy/Oa, 
Oz/Op1 02/OpNR. 02/004 
则 
Pi1=TPxT. (19) 
和 矩阵 I 中 的 后 3 行 可 通过 动力 学 映射 得 到 . 显 
然 , (16) 式 成 立 的 条 件 是 Y 一 0, 这 需要 较 好 的 z 初 
始 猜测 值 , 但 这 在 一 般 情 形 下 不 能 简单 获得 . 在 本 


文 工 作 中 , 我 们 借助 了 上 其 有 全 局 优化 性 质 的 粒子 群 


算法 解决 此 问 
两 步 , 即 首先 用 粒子 群 
化 解 , 然后 再 针对 此 解 
局 部 最 优 解 . 


3.2 


假设 


(ad, 63, 49, 63), 


假定 已 知 第 1 个 观测 
度 pi 和 1, 在 本 问题 中 种 


X = 


待 优化 的 成 本 函数 仍 | 


题 . 整个 双 r 算 法 初 轨 确 定 的 流程 分 


算法 求解 (12) 式 寻找 全 局 优 
采用 (16) 式 逐步 夫 代 以 获得 


本 文 所 用 可 行 域 方法 的 改进 
LAS NK (> 2) 组 TSA 资 料 


d=1,2,---,NK. 


资料 的 站 心 距离 及 其 速 
村 估量 即 为 


(p1, p1)". 


日 可 以 表示 成 (12) 式 的 形式 ， 


但 每 个 时 刻 的 观测 残 差 变 为 


Ya = [cos ôg (Qa 一 


aa) (64 一 


ba): 


cos63(63 — åa) (dg — 09)", d > 2. 


整个 求解 过 程 与 上 述 双 7 方法 类 似 , 即 首先 利 


j 智 能 优化 算法 在 可 行 


ChinaXiv 合 作 期 刊 


—p1 Sin 01 cos Q1 0 0 
— pı sind; sina, 0 0 
pı COS 01 0 0 
0% /06, DZ/DaNR DZ/OONR 
Oy /061 Oy /Oanr Oy /OOnR 
02/00, 0z/Oanr DZ2/OONR 
FMB AR nC. 为 计算 待 估量 的 协 方差 
ERE, 引入 协 方差 矩阵 
Py, 0 0 
0 P ee 0 
Py =EYY"]=| | “i l . , 
0 0 Py... 1 


WU 4 (rH at AN Bb Fr 2 FB ME 5 IY AODA. E 
样 ， 对 每 个 时 刻 的 协 方 差 矩阵 Py 而 言 除 需 考 
虑 该 时 刻 该 观测 量 自身 的 方差 ， 还 需 考虑 来 
自 第 1 个 TSA 的 观测 量 的 误差 贡献 ， 即 形式 仍旧 
如 (18) 式 所 示 , 其 中 


cos” dga2, 0 0 0 
po _ 0 oS, 0 0 
a 0 0 cos?59a2 0 
0 0 0 a 
为 观测 量 自身 的 方差 信息 ， 
o2 0 0 0 
(2) _ dYa 0 of 0 0 
Ya apd | 0 dB 0 
0 0 0 oF 
(az, 62, a2, 07) 
为 第 1 个 TSA 的 方差 信息 映射 至 后 续 观 测 时 刻 的 方 


区 域内 全 局 寻 优 ， IRE 


的 怀 使 得 残 差 (12) 式 最 小 ， 其 次 在 局 部 区 域 采 
(16) 式 进一步 迭代 求解 . 与 双 r 方 法 的 不 同 之 处 在 


ik 


差 贡 献 . EA AHERE 


dYa i 
(aF, 69, a9, 67) 
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可 通过 初始 时 刻 动力 学 映射 简单 求解 . 

与 双 7 方 法 类 似 , (17) 式 仅 给 出 了 待 估量 的 协 方 
ZEEE, AE RN AVG 速度 的 协 方差 信息 ， 
需 进一步 求解 ， 计 算 公 式 仍 如 (19) 式 所 示 , (AX AY 
定义 变 为 


x = (X, a9, ô, ation) = (P1, 61,08, 69, a9, 62)". 


矩阵 已 总 的 定义 为 
Py 0 0 0 0 
0 o, 0 0 0 
Py= 0 co 0 0 
0 


0 
0 0 0 æ 
0 


过 (11) 式 几何 关 


此 (19) 式 中 的 Jaccobi 和 矩阵 工 可 通 


3.3 ”改进 Laplace 算 法 的 协 方差 矩阵 计算 


ENR (> 3) 组 观测 量 , 引入 各 观测 时 刻 的 “ 残 
xB” 


(F* xo + G% to — Xm)Um— 
(Fizo + Gmžo — 2m) Am 
(Fyo + Ginto — Ym)Ym— 
(Fizo + Gmo — Zm) Um 
(FX £o + G>, to — Xm)Um— 
(FnYyo + Gm to — Ym)Am 


Ym = 


此 方法 中 待 估量 应 修改 为 


所 有 的 残 差 向 量 为 


Yı 


E RRRA AGN YT TIE A Wr TIE, EXB 
以 及 工 矩 阵 进行 了 相应 变化 . 定义 
Oy,, Ol Anes Hm, Vm) 


Ym = 


Ol Anns Um, Vm) Alam, m) ` 
工 形 如 
Yı 
T= 
YNR 
不 难得 到 
a, 0 0 0 
0 of 0 0 
Py =T : : rt 
0 0 Gene O 
0 0 se 0 Sy 


(20) 

代入 (17) 式 即 可 计算 待 估量 的 协 方差 矩阵 . 
3.4 ”粒子 群 算法 

在 本 文 所 作 工 作 中 , 对 于 每 个 方法 所 求 参数 的 
寻找 ,比如 双 7 方 法 中 的 p1 和 pnwr, 可 行 域 方法 中 的 
Pp1 和 ,是 采用 粒子 群 算法 进行 直接 搜寻 , 使 得 目 
这 里 简要 介绍 粒子 群 算法 以 及 其 计算 

eh 
ee sl on 
主 自己 距离 栖息 地 最 近 的 位 置 并 reels 
信息 . 即 每 个 粒子 根据 自己 的 位 置 和 速度 进行 移动 ， 
根据 最 优 值 的 信息 进行 速度 的 修正 . 其 中 对 速度 的 
修正 考虑 权重 、 加 速 因子 的 修正 公式 . 具体 关于 粒 
子 群 优化 算法 的 计算 过 程 如 下 4. 
WHR: 假设 有 / 维 需要 优化 的 变量 ， 每 一 个 
维度 的 大 小 范围 为 [au bal, 设 定 的 初始 种 群 数 size- 
pop， 即 粒子 数 . 保证 每 个 维度 在 设 定 的 数值 范围 
内 , 生成 sizepop 个 初始 粒子 , 每 一 个 粒子 即 可 计算 
得 到 初始 适应 度 fitness, 在 本 问题 中 即 为 观测 残 差 ， 
另外 根据 飞行 最 大 速度 限制 , 同时 生成 该 粒子 对 应 
的 速度 ; 

步骤 2: 第 ! 次 迭代 过 程 , 即 第 ! 代 种 群 , 根据 公式 
依次 对 权重 、 加 速 因子 、 速 度 以 及 粒子 更 新 , 需要 
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考虑 到 的 是 速度 更 新 后 一 旦 超过 速度 限制 范围 , 偏 
大 即 取 速 度 最 大 值 ， 偏 小 即 取 速度 最 小 值 . 而 一 旦 
粒子 更 新 后 超出 初始 设 定 范围 , 重新 在 粒子 约束 范 
围 内 寻找 一 个 随机 数 作为 其 值 , 并 且 对 于 该 粒子 速 
度 更 新 为 0, 重新 计算 适应 度 . 再 进行 个 体 极 值 和 群 
Be 该 次 迭代 计算 完成 ; 

又 3: 记录 每 一 代 最 优 取 值 以 及 粒子 位 置 记 
录 . a 
数 ; 

步 又 4: 根据 记录 的 每 一 代 最 优 值 , 找到 最 优 值 
以 及 对 应 的 粒子 , 求解 结束 . 

对 于 本 问题 的 计算 , 即 优 化 变量 为 需要 猜测 的 
值 如 mr 和 pm pf, 适应 度 的 计算 即 目标 函数 的 
计算 , 而 优化 变量 的 范围 、 种 群 数 和 迭代 次 数 的 选 
取 则 需要 根据 经 验 和 实际 情况 给 出 . 


4 数据 处 理 与 试验 结果 
4.1 数据 准备 
本 文中 采用 数据 参考 国际 天 文联 合 会 (Inter- 
national Astronomical Union, IAU) 所 属 的 小 天 体 
中 心 (Minor Planet Center，MPC) 网 站 的 真实 观 
测 数据 以 及 测 站 数据 . 利用 ESA (European Space 
Agency) 维 护 的 NEODyS (Near Earth Objects Dy- 
namic Site) 所 给 的 小 天 体 Toutatis (4179) 的 轨道 信 
息 作 为 积分 初 值 预 报 短 弧 段 , 给 出 该 弧 段 内 观测 时 
刻 对 应 的 赤 经 和 赤 纬 及 其 变 率 , 并 赋 以 随机 误差 ， 
赤 经 和 赤 纬 的 误差 不 超过 1”， 其 变 率 的 误差 不 起 
过 1(”) - h7. 
4.1.1 Toutatis (04179) 基 本 信息 
表 1 给 出 了 Toutatis (04179) 的 半 长 径 、 偏 心 
率 、 对 应 简约 儒 略 日 以 及 周期 的 基本 信息 . 
4.1.2 角度 变 率 误差 估计 
以 赤 经 为 例 , 假设 采用 线性 方法 给 出 观测 角度 
的 变 率 
= ʻa z änt = — 3 Š läst. (21) 
取 一 阶 项 近似 计算 , 则 (最 后 一 项 为 系统 差 ) 
ÔQ — O01 Aa 
At (At? 


5(At) 4 sat, (22) 


6a = 


*http://www.iausofa.org/ 
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因此 协 方差 为 

B[oada] = m Ani (At) ' Laat 
=a [AH Laat 


假定 各 观测 时 刻 都 很 准确 ( 即 上 式 中 的 5(Ab) 
是 一 小 量 , 在 目前 的 时 间 单位 系统 下 这 是 一 个 可 接 


受 的 假设 ), 且 假 定 E(da2) = E(da,) = 0. 则 上 式 
E[6ada] ~ ra lB (Ganson) aodan 
7a (At)? . (24) 


在 真实 计算 中 , HAM ARA REAL, 为 确保 
赤 经 和 赤 纬 变化 速率 的 准确 性 , 假设 两 次 观测 之 间 
时 间 差 At 约 1.4 h. WR a ee 
个 小 量 , 对 于 赤 经 约 为 4(”) . h-2， 与 第 1 项 量 级 相 
比 可 忽略 不 计 ， 则 此 时 可 估算 到 赤 经 以 及 赤 纬 变 率 
残 差 约 为 1(”) hol. 在 本 文 的 工作 中 将 以 该 值 作为 
随机 残 差 加 入 . 


表 1 Toutatis (04179) 基 本 信息 . MJD 表 示 简 约 儒 略 日 ， 


a 表 示 半 长 径 , e 表 示 偏 心率 , 了 表示 周期 . 
Table 1 Basic information of Toutatis (04179). 


MJD represents Modified Julian Date, a 


represents the semi-major axis, e represents the 


eccentricity, and T represents the period. 
MJD a/AU e T/d 


59800.0 2.54 0.62 


1482.68 


4.1.3 测 站 速度 计算 

由 于 观测 时 的 坐标 系 为 地 心 天 球 坐 标 系 ， 而 
测 站 位 置 为 地 固 坐 标 系 下 的 位 置 , 所 以 会 考虑 地 
心 天 球 坐 标 系 与 地 固 坐 标 系 之 间 的 转换 ， 本 文采 
用 Sofa! 提 供 的 程序 进行 坐标 转换 . 另外 , 在 计算 赤 
经 和 赤 纬 变 率 时 , 需要 计算 测 站 的 速度 信息 , 而 测 
站 位 置信 息 表达 式 为 


= (HG)T Yr, 5 
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其 中 ,及 表示 在 地 心 历 元 平 赤道 下 的 测 站 位 置 矢 
量 , 表示 在 地 固 系 下 的 测 站 位 置 矢 量 , (HG)T 表 示 
上 述 两 坐标 系 之 间 的 转换 矩阵 ， 即 包括 地 球 自转 、 
极 移 、 岁 差 和 章 动 . 若 需要 求 测 站 速度 , 则 需要 对 
上 和 式 进行 求 导 , 而 (HG) 德 阵 的 导数 由 理论 公式 导 
出 较 复 杂 . 故 采 用 数值 方法 一 Lagrange 四 阶 插值 公 
式 直 接 对 计算 的 测 角 资 料 进行 插值 , 并 对 该 插值 公 
式 进 行 求 导 , 给 出 测 角 资料 的 变 率 . Lagrange 四 阶 
MEZA RAN F: 
_ (æ—z)(z 
FO) = Ee 
(x — zo) (x — zx2)(7 — z3) 
(£1 — £o) (£1 — £2) (£1 — T3) 
(x — zo) (x — z2) (x — z3) 
(£1 — £o) (£1 — £2) (£1 — T3) 
(x — zo) (x — xı) (z — z2) 
(£3 — £o) (£3 — x1)(£3 一 va l (25) 
考虑 观测 时 刻 为 tf, 时 间 差 At, Wao. zis wey 
£33) A At—2At, t—At. t+At, t+2At, 对 (25) 式 
进行 求 导 , 可 得 
f(z) = 
41.4 光 行 差 计 算 
在 进行 轨道 积分 时 ， 我 们 采用 了 RKF78 积 分 


器 , 由 于 时 间 长 度 很 短 , 采用 力 模 型 只 考虑 太阳 的 
引力 . 对 于 人 造 卫 星 而 言 , 由 于 与 测 站 距离 较 短 , 初 


ar 


X2)(x — x3) 


Yo 


8y3 — 8y2 — Yı + Yo 


12At 


(26) 


Tk true 
a T y O(10-2) . 


C 


该 时 间 差 造成 的 角度 误差 约 0.0001”, 可 忽略 不 计 . 
4.1.5 其 他 说 明 

粒子 群 算法 在 生成 种 群 时 , 需要 对 优化 变量 进 
行 限 制 , 并 设置 迭代 次 数 和 种 群 数 , 对 于 可 行 域 方 
法 , 文中 设置 种 群 数 为 1000 个 ， 迭代 次 数 为 300 次 ， 
视 向 距离 变化 范围 为 (ap,2AU), 视 向 距离 的 变 率 
变化 范围 为 (一 0.00593 AU - s7}, 0.00593 AU - s71); 
而 对 于 双 7 方 法 ,设置 种 群 数 为 1000 个 ， 迭 代 次 数 
为 300 次 , 视 向 距离 的 范围 为 (as,2 AU). 

本 文 所 用 弧 段 数据 请 见 附录 , 为 与 真实 轨道 值 
进行 对 比 ， 本 文保 留 真实 观测 数据 对 应 的 时 刻 以 
及 测 站 信息 ,而 其 观测 信息 ( 赤 经 和 赤 纬 以 及 其 变 
率 ) 由 拟定 轨道 进行 轨道 积分 再 加 上 随机 误差 给 出 . 
4.2 ”试验 结 

在 本 部 4 eee. 双 r 方 法 、 
可 行 域 方法 共 3 种 方法 在 不 同 弧 段 情形 下 的 计算 结 
果 . 在 本 文中 , 将 少 于 一 天 观测 时 长 、 观 测 数目 少 
于 5 的 弧 段 定义 为 极 短 弧 , 超过 150 d 的 弧 段 定义 为 
KIE, 其 余 弧 段 定义 为 普通 弧 段 . 在 本 工作 中 , 为 
了 防止 随机 误差 偶然 影响 对 结果 的 判断 , 故 对 每 一 
组 原始 弧 段 数据 多 次 加 入 新 的 随机 误差 , 生成 不 同 
的 弧 段 数据 再 进行 计算 , 并 比较 最 后 结果 . 另外 , 为 
了 防止 由 于 力学 模型 造成 的 结果 差异 影响 其 结果 


轨 确 定 可 忽略 光 行 差 ; 但 是 对 于 近 地 小 天 体 而 言 ， 
与 测 站 距离 较 远 , 甚至 可 以 达到 1 AU, 根据 观测 资 
料 t;, 此 时 小 天 体 与 测 站 距离 p;, 时 间 修 正 为 

假设 pj; 为 1 AU. bt; ~ 500 s， 如 果 小 天 体 运 
动 速度 为 20 km:s-!， 此 时 造成 观测 角度 误差 约 为 
10-5 rad, 相当 于 2”, 因此 需要 进行 光 行 差 修正 . 实 
际 上 , 时 间 修 正 公 式 应 该 为 


i P , 

其 中 c 表 示 光 速 . 上 式 应 该 迭代 求解 , 但 是 

[ri(ts — 867") — ri (ti) 
Cc 


bt; — ot" = 


与 双 7 方 法 和 改进 Laplace 方 法 的 对 比 , 可 行 域 方法 

在 二 体 问题 下 进行 积分 

4.2.1 第 2 步 迭 代 的 讨论 
在 本 文 工 作 中 ， 当 采 


J 可行 域 方法 或 双 r 方 


法 时 ， 理 论 上 求解 过 程 需 分 两 部 完成 第 1 步 
是 采用 智能 优化 算法 尽 可 能 寻找 全 局 的 优化 解 ， 


第 2 步 是 在 第 1 步 的 基础 上 进一步 采用 迭代 方式 ( 参 
见 (16) 式 ) 以 期 望 对 全 局 优化 解 进行 进一步 优化 ( 即 
使 得 观测 残 差 进 一 步 变 小 ). 然而 我 们 实际 的 数值 
验证 表明 , 很 多 时 候 采 用 第 2 步 的 局 部 迭代 反而 会 
使 得 残 差 增 大 . 造成 这 种 现象 的 原因 是 在 第 2 步 迭 
代 的 公式 推导 中 假设 Y 一 0, 去 掉 了 二 阶 导 数 项 
YTƏB/ƏX, 这 在 精密 定 轨 中 是 常见 的 . 但 在 极 短 
弧 定 轨 中 ， 由 于 缺乏 足够 的 信息 约束 轨道 ,因此 
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Y 一 0 的 假设 可 能 不 再 成 立 , 也 即 (15) 式 得 到 (16) 
式 过 程 中 的 假设 


不 再 成 立 . 基于 这 种 现象 , 我 们 在 极 短 弧 定 轨 的 情 
形 下 实际 上 仅 进 行 了 第 1 步 并 将 其 结果 作为 最 终 定 
轨 结 果 输 出 . 
4.2.2 ” 极 短 弧 段 试验 结果 
图 4 和 图 5 中 展示 了 极 短 弧 段 的 计算 结果 , 该 弧 
段 时 长 约 2 h,， 共 4 条 观测 数据 . 这 里 生成 了 20 组 测 
试 数据 (每 组 测试 数据 含 4 条 观测 值 , 但 每 组 数据 中 
的 观测 值 施 加 了 不 同 的 随机 差 ) 横 坐 标 代表 测试 
数据 的 组 数 , 纵 坐标 代表 每 一 组 测试 数据 计算 半 长 
径 与 真实 半 长 径 之 差 或 残 差 RMS. 在 图 中 未 标示 出 
利用 双 r 方 法 计算 无 结果 的 情况 . 在 20 组 试验 中 , 改 
进 Laplace 方 法 都 失败 ， 数 值 验算 时 发 现 会 出 现 迭 
代 解 在 两 个 点 之 间 来 回 跳 转 的 情况 , 利用 高 斯 方法 
直接 求解 八 阶 方程 可 以 证 实 这 种 现象 由 上 述 多 解 


问题 造成 . 


3 
25| a | 
22 Admissible Region Method 
at < Double-r Method 
=; 
< 
15 t f 
S >< 
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~% Y * 
0 * L = 1 
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Group Number 
A4 弧 段 1 半 长 径 计算 结果 


Fig.4 The semi-major axis result of arc 1 


在 极 短 弧 结合 可 行 域 方 法 的 情形 下 ， 本 身 弧 
段 的 几何 状况 可 能 会 导致 寻 优 之 后 给 出 的 轨道 出 
现 半 长 径 极 大 或 者 偏心 率 趋 于 1 的 情况 ， 此 时 该 轨 
道 仍 在 2.3 节 所 描述 的 可 行 域内 , 但 显然 不 符合 实 
际 情况 .在 本 文 的 工作 中 ,我 们 人 为 添加 了 半 长 
径 a < 5.2 AU 的 约束 . 即使 如 此 ,可 行 域 方 法 的 定 
轨 结 果 中 有 多 组 算 例 出 现 半 长 径 o 接 近 5.2AU 的 情 
%, 但 其 残 差 ( 见 图 中 仍旧 在 1 以 内 . 在 极 短 弧 段 时 ， 
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定 初 轨 的 计算 结果 很 大 一 部 分 取决 于 轨道 几何 约 
R, 需要 加 入 更 多 的 约束 才能 得 到 更 合理 的 结果 . 
总 的 来 说 , 对 于 极 短 弧 段 , 改进 Laplace 方 法 基本 不 
适用 ， 而 可 行 域 方法 和 双 r 方 法 虽然 大 多 数 都 能 给 
出 定 轨 结 果 , 但 由 图 4 可 知 这 些 结果 可 信和 度 较 低 . 若 
需要 计算 获得 更 准确 的 轨道 , 需要 加 入 更 多 的 约束 . 


0.3 
0.25| > “x Admissible Region Method 
1 X Double-r Method 
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米 
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0.1 E a Pe 
0.05 
o =Z HIATT RI NA 
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图 5 弧 段 1 残 差 计算 结果 
Fig.5 The RMS result of arc 1 
4.2.3 普通 长 度 弧 段 试验 结果 


图 6 和 图 7 中 , 考虑 弧 段 2, 时 长 约 5 d, 观测 数目 
为 4 个 , 与 弧 段 1 相 比 只 是 弧 段 时 长 增加 , 生成 多 组 
数据 比较 计算 得 到 的 半 长 径 与 残 差 . 


2.5 
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Ae 弧 段 2 半 长 径 计 算 结果 


Fig.6 The semi-major axis result of arc 2 
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Z Admissible Region Method 改善 , 但 双 7 方 法 的 残 差 结果 反而 变 差 , 半 长 径 结果 
ate ita Metod 与 弧 段 2 近似 . 双 r 方 法 产生 该 结果 的 原因 将 在 后 续 
ta 进行 分 析 . 
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图 7 弧 段 2 残 差 计 算 结 果 ne 
X y * 六 % $ 
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图 7 中 可 行 域 方 法 的 残 差 明显 高 于 其 他 两 种 方 
法 , 但 大 部 分 仍旧 小 于 1”, 而 图 6 中 由 改进 Laplace 
方法 和 双 7 方 法 给 出 的 计算 结果 与 真实 半 长 径 非 常 
吻合 , 并 且 这 两 种 方法 计算 所 得 的 结果 非常 相似 . 
该 弧 段 结果 与 弧 段 1 进行 对 比 , 在 观测 弧 段 时 长 增 
长 的 情况 下 , 改进 Laplace 方 法 和 双 r 方 法 能 给 出 较 
好 的 结果 , 可 见 利用 这 两 种 方法 在 处 理 初 轨 计 算 时 ， 
弧 段 不 能 太 短 ; 而 可 行 域 方法 的 结果 与 弧 段 1 近似 ， 
其 计算 结果 与 弧 长 无 明显 关系 . 对 于 可 行 域 方 法 观 
测 数目 相同 、 弧 长 变 长 , 但 结果 却 不 尽 人 意 的 情况 ， 
猜测 可 能 是 有 两 种 原因 造成 的 , 一 是 由 于 赤 经 和 赤 
“ABR (ae, 上) 误差 本 身 较 大 ,二 是 约束 不 够 充分 ， 
可 行 域 方 法 受到 观测 数目 的 影响 更 大 而 非 观 测 弧 
长 . 为 验证 这 两 种 猜想 , 进行 了 两 个 试验 , 一 是 在 进 
行 粒子 群 优化 算法 全 局 搜索 时 , 不 对 赤 经 和 赤 纬 变 
率 进 行 约束 , 仅仅 将 第 1 组 赤 经 和 赤 纬 变 率 作为 输 
A, 试验 发 现 结果 并 没有 明显 变化 , 故 排除 第 1 种 猜 
想 ; 二 是 增加 观测 数目 , 仍旧 保持 弧 长 约 5 d, 观测 
数目 变 为 15, 即 弧 段 3 进行 计算 . 图 8 和 9 中 分 别 展示 
了 弧 段 3 半 长 径 和 残 差 的 计算 结果 . 对 比 图 6 和 7 的 
结果 不 难 发 现 , 此 时 可 行 域 方法 结果 明显 有 所 改善 
另 一 方面 双 r 方 法 出 现 半 长 径 与 真实 值 相 近 但 残 差 
结果 较 大 的 现象 , 改进 Laplace 方 法 在 3 者 中 仍然 表 
现 最 好 . 

其 次 ， 考虑 时 长 约 6 d， 观 测 数目 为 22 个 的 弧 
段 4. 图 10 和 11 中 给 出 弧 段 4 的 多 组 计算 结果 . 对 于 
弧 段 4, 改进 Laplace 方 法 和 可 行 域 方法 都 能 给 出 较 
好 的 结果 . 相 比 于 弧 段 2, 可 行 域 方法 计算 结果 有 所 
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8 ” 弧 段 3 半 长 径 计算 结果 


Fig.8 The semi-major axis result of arc 3 
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图 9 ” 弧 段 3 残 差 计算 结果 


Fig.9 The RMS result of arc 3 


最 后 , 我 们 考虑 了 弧 段 5, 时 长 21 d, 共 62 个 观 


测 数据 点 . 图 12 和 图 13 中 ,针对 弧 段 5 生成 多 组 仿 
真 观测 数据 , 并 将 定 轨 给 出 的 半 长 径 和 残 差 与 之 前 


弧 段 的 结果 进行 
法 的 计算 结果 非 


比较 . 如 图 13, 首先 改进 Laplace 方 


常 稳定 , 半 长 径 基 本 在 一 个 小 区 域 


范围 内 波动 , 说 明 在 观测 数据 更 多 时 约束 更 多 , 计 
算 结果 能 收敛 到 固定 区 域 , 但 此 时 半 长 径 与 真实 


的 1988 年 12 月 26 


结果 之 差 (图 12) 相 较 于 前 几 条 弧 段 更 大 . 后 经 查证 ， 
该 组 数据 的 时 间 起 始 于 1989 年 1 月 3 日 ， 而 在 8 d 前 


H, 该 天 体 与 地 球 发 生 了 近 距 离 的 


交会 . 选取 1989 年 1 月 4 日 开始 的 62 个 观测 样本 点 ， 
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总 时 长 约 19.6 d, RAIG, 进行 初 轨 计 算 结 果 如 真实 值 相差 较 小 , 并 且 残 差 保 持 在 才 以 下 . 
图 14 和 15 所 示 ， 此 时 改进 Laplace 方 法 的 半 长 径 与 
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Fig.12 The semi-major axis result of arc 5 Fig.13 The RMS result of arc 5 
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另外 , 针对 于 上 述 的 弧 段 1-5, 随 着 弧 段 时 长 和 
观测 数目 的 增加 ,可 以 看 到 3 种 方法 计算 得 到 的 半 
长 径 与 真实 半 长 径 之 差 越 来 越 小 , 可 见 观测 约束 越 
多 , 得 到 的 轨道 信息 越 真实 . 然而 此 时 残 差 结 果 并 
未 表现 出 这 种 现 P 甚至 某 些 弧 段 和 方法 的 残 差 计 
算 结 果 还 会 增 大 , 可 见 在 初 轨 计 算 中 , 残 差 结果 只 
是 参考 项 , 并 未 起 决定 性 作用 . 

针对 弧 段 4、5、6 对 于 双 7 方 法 的 计算 结果 , 其 
残 差 不 理 想 , 但 是 其 半 长 径 与 真实 半 长 径 相 近 , 这 
里 给 出 双方 法 关于 半 长 径 和 残 差 收敛 区 域 示 意图 
如 图 16 和 17. 图 16 描 述 了 计算 半 长 径 与 真实 半 长 径 
之 差 的 函数 ba(p1, pr), 单位 为 AU, 图 17 描 述 了 计算 
残 差 的 函数 RMS(p1, pz?)， 单 位 为 ", 本 问题 即 在 上 
述 两 图 中 寻找 同时 满足 极 小 的 点 . 如 图 17, 尽管 展 
示 了 残 差 小 于 80/ 的 区 域 , 但 其 收敛 域 仍然 只 存在 
于 极 细 的 直线 上 ; 反观 图 16, 关于 半 长 径 极 宽 的 收 
DX SR. 当 目 标 函 数 为 残 差 时 , 很 难 找到 全 局 
极 小 点 ， 导 致 得 到 的 残 差 较 大 , 但 其 位 置 所 对 应 的 


半 长 径 却 与 真实 值 接近 
0 1 
0.05 poe 
= 0.6% 
< 0.1 = 
a 3 
Q 0.4<1 
0.15 52 
0.2 0 
0.1 0.15 0.2 
pı/ AU 


|16 ” 双 7 方 法 关于 半 长 径 收 敛 区 域 示意 图 


Fig.16 The converging region about the semi-major axis of 
double-r method 


5 讨论 

当 弧 上 段 较 长 时 , 本 文 基于 时 间 级 数 展开 方式 的 
改进 Laplace 方 法 适用 性 会 受到 挑战 ， 因 为 此 时 该 
级 数 有 可 能 不 收敛 . 但 双 r 方 法 仍然 适用 . 我 们 开 
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展 了 几 组 数值 验证 , 结果 如 表 2 所 示 . 可 以 看 到 对 于 
弧 段 7 和 弧 段 8, 计算 出 的 轨道 与 真实 轨道 根 数 非常 
相似 , 但 是 弧 段 9 的 结果 比 前 两 个 弧 段 较 差 (可 能 是 
此 时 时 间 过 于 长 , 三 体 问题 下 的 轨道 积分 和 NN 体 问 
题 (N-body Problem) 下 的 轨道 积分 有 了 较 大 差异 ). 


0:12 80 


0 
0.1 0.105 0.11 0.12 


pı/ AU 


图 17 Wr ATER RAIS KANE E 


0.115 


Fig.17 The converging region about the RMS of double-r 
method 


一 般 来 说 , 根据 定义 的 马 氏 距离 可 进行 轨道 关 
联 的 计算 , 相当 于 是 “ 聚 类 ”的 问题 . Milani $2 
到 了 可 行 域 方法 在 轨道 关联 中 的 应 用 ，Gronchi 
等 25 利 用 轨道 积分 进行 初 轨 确定 的 同时 也 进行 了 
轨道 关联 的 工作 ， 也 是 利用 了 可 行 域 方法 ， 除 此 
之 外 还 有 许多 学 者 依据 可 行 域 的 思想 针对 绕 地 目 
标 进行 了 轨道 关联 的 工作 . Ti Gooding! ATH HAY 
双 7 方 法 的 思想 , 其 实 也 可 以 在 轨道 关联 中 应 用 , 即 
将 一 些 相隔 时 间 较 长 的 观测 资料 (数目 大 于 2) 放 在 
一 起 进行 初 轨 计 算 , BRAK BE ANAL F A 
一 观测 目标 的 弧 段 . 


6 结果 

本 文 针 对 偏心 率 较 大 的 天 体 Toutatis, 选取 少 
于 一 天 观测 时 长 的 极 短 弧 段 、 普 通 长 度 弧 段 以 及 
超过 150 d 时 长 的 长 弧 段 ,分 别 采用 3 种 方法 (改进 
Laplace 方 法 、 双 7 方法 和 可 行 域 方 法 ) 进 行 初 轨 计 
算 . 
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表 2 双 7 方 法 应 用 于 长 弧 段 的 计算 结果 . a、e、it、Q、w 和 MM 分 别 表示 轨道 6 要 素 中 半 长 径 、 偏 心率 、 轨 道 倾角 、 升 交点 经 
度 、 近 心 点 幅 角 和 平 近 点 角 . 


Table 2 The results of the long arcs using double-r method. a, e, i1, Q, w, and M respectively represent 


the six orbital elements: semi-major axis, eccentricity, inclination, longitude of the ascending node, 


argument of perigee, and mean anomaly. 


Arc i 8 9 
Arc length/d 149.58 199.16 250.21 
Number of observation data 91 91 289 
Aa/AU 1.8x1073 —5.1x1073 —5.1 x 107% 
Ae 2x 1074 —0.7 x 107? 一 5.7 x 1078 
Ais /° 5.5 x 107% 3.6 x 107% —5.7 x 107% 
AQ/° —7.3 x 107? -1x10-* 5.7x 10 
Aw/° —5.1 x 107? 1.2 x 107? —0.2 
AM/° —3.01 x 107? —9.9 x 107? 2.9 


结果 表明 , 对 于 极 短 弧 段 , 尤其 是 观测 资料 少 
的 情况 , 可 行 域 方 法 和 双 7r 方 法 都 能 给 出 残 差 较 小 
的 结果 ,而 改进 Laplace 方 法 通常 无 法 给 出 计算 结 
R, 如 若 需 要 给 出 符合 真实 情况 的 轨道 , 需要 加 入 


更 多 约束 进一步 确定 轨道 ; 对 于 普通 的 弧 段 , 我 们 
的 研究 表明 双 r 方 法 和 改进 Laplace 方 法 的 计算 结 
果 相对 更 依赖 于 观测 弧 段 长 度 而 非 观 测 数目 , 而 可 
行 域 方 法 结果 在 观测 弧 段 相同 时 , 更 依赖 于 观测 数 
目 ; 在 观测 数目 较 多 时 , 改进 Laplace 方 法 的 结果 精 
度 和 计算 速度 更 突出 ; 对 于 较 长 的 弧 段 , 双 r 方 法 具 
备 可 行 性 , 可 在 近 地 天 体 观测 资料 的 轨道 关联 中 有 
所 应 用 . 

在 针对 不 同 弧 段 时 长 和 观测 数目 的 情况 下 , 可 
以 采取 不 同 的 方法 进行 初 轨 确 定 , 在 弧 段 极 短 时 选 
取 双 7 方法 或 者 可 行 域 方法 ; 对 于 普通 长 度 至 几 天 
长 度 的 弧 段 , 选取 改进 Laplace 方 法 , 此 时 如 若 选 择 
其 他 方案 , 观测 数目 较 多 选择 可 行 域 方法 , 观测 数 
目 较 少 选择 双 r 方 法 ; 对 于 长 弧 段 可 以 选择 双 7r 方 法 . 


I 
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Initial Orbit Determination Based on Intelligent Optimization 
Algorithm 


LIU Xinl23 HOU Xiryunl23 LIU Lint?3 GAN Qing-bo* YANG Zhi-tao* 
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Apstracr Classical methods for initial orbit determination (IOD) include Laplace method, Gauss 
method and their variations. In addition to this, based on the characteristic of optical observation data 
nowadays, experts propose some other IOD methods, like Double-r method and admissible region method. 
One of the ways to determinate the orbit through double-r method is to guess distances of the target from 
the observer at two epochs—usually at the first and the last one. By doing so, we can solve the Lambert 
problem and use its solution as the initial guess of the orbit. Furthermore, we can improve the initial guess 
by iterations to reduce the root mean square (RMS) of the observations. The admissible method is based 
on the concept of attributable (longitude, latitude and their rates). With some conceptions, the admissible 
region described by the range and range rate from the observer is characterized. Using triangulation we 
can find the nodal point that makes the RMS minimal. In our work, we apply one intelligent optimization 
method—the particle swarm optimization method to the two methods, based on simulated and real data, 
and compare the results with that of modified Laplace method. At last, we briefly discuss the possibility 
of applying the double-r method to the orbit link problem. 


Key words celestial mechanics, minor planets, methods: numerical 
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附录 


表 3 各 弧 段 初 值 


‘Table 3 The initial values of arcs 


Arc Number af? 6/° afl) <a> BO .s-!) Date (TDB*) 
1 54.888 17.836 0.12133 0.03258 1989-01-04 19:30:56.18 
2 55.010 17.868 0.12097 0.03216 1989-01-04 20:31:05.97 
! 3 55.056 17.881 0.12090 0.03200 1989-01-04 20:53:56.28 
4 55.176 17.912 0.12085 0.03159 1989-01-04 21:53:56.56 
1 299.250 —20.476 —0.01549 —0.00293 1988-07-12 08:17:26.03 
2 298.925 —20.540 —0.01577 —0.00310 1988-07-13 05:57:40.91 
i 3 297.737 —20.768 —0.01718 —0.00307 1988-07-16 08:43:55.79 
4 297.341  —20.843 —0.01757 —0.00311 1988-07-17 08:28:26.13 
1 49.880 16.466 0.13272 0.03920 1989-01-03 05:12:11.54 
2 49.895 16.470 0.13272 0.03914 1989-01-03 05:19:11.45 
3 54.889 17.836 0.12134 0.03259 1989-01-04 19:30:56.18 
4 55.011 17.869 0.12098 0.03216 1989-01-04 20:31:05.98 
5 55.057 17.881 0.12090 0.03200 1989-01-04 20:53:56.28 
6 55.177 17.913 0.12085 0.03159 1989-01-04 21:53:56.57 
7 58.566 18.754 0.11190 0.02726 1989-01-06 02:14:55.70 
3 8 58.964 18.849 0.11252 0.02601 1989-01-06 05:48:08.09 
9 60.786 19.268 0.10635 0.02415 1989-01-06 22:04:55.80 
10 60.813 19.274 0.10636 0.02406 1989-01-06 22:19:56.09 
11 61.300 19.383 0.10476 0.02446 1989-01-07 02:42:56.18 
12 61.318 19.387 0.10464 0.02439 1989-01-07 02:52:55.80 
13 61.608 19.453 0.10353 0.02324 1989-01-07 05:40:25.85 
14 61.625 19.457 0.10352 0.02317 1989-01-07 05:50:26.33 
15 64.878 20.139 0.09491 0.01917 1989-01-08 13:37:25.85 
1 49.879 16.466 0.13272 0.03920 1989-01-03 05:12:11.54 
2 49.895 16.470 0.13272 0.03914 1989-01-03 05:19:11.45 
3 54.889 17.836 0.12134 0.03259 1989-01-04 19:30:56.18 
4 4 55.010 17.869 0.12098 0.03216 1989-01-04 20:31:05.98 
5 55.056 17.881 0.12090 0.03200 1989-01-04 20:53:56.28 
6 55.177 17.913 0.12085 0.03160 1989-01-04 21:53:56.57 
7 58.566 18.754 0.11190 0.02726 1989-01-06 02:14:55.70 


*Barycentric Dynamical Time 
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表 3 续 
Table 3 Continued 
Arc Number a/® 6/° àa :sl) 6/((") 7") Date (TDB) 
8 58.964 18.849 0.11252 0.02601 1989-01-06 05:48:08.09 
9 60.786 19.268 0.10635 0.02415 1989-01-06 22:04:55.80 
10 60.812 19.274 0.10636 0.02406 1989-01-06 22:19:56.09 
11 61.300 19.383 0.10476 0.02446 1989-01-07 02:42:56.18 
12 61.317 19.387 0.10464 0.02439 1989-01-07 02:52:55.80 
13 61.608 19.453 0.10353 0.02324 1989-01-07 05:40:25.85 
14 61.625 19.457 0.10352 0.02318 1989-01-07 05:50:26.33 
4 15 64.878 20.139 0.09491 0.01917 1989-01-08 13:37:25.85 
16 64.888 20.141 0.09492 0.01913 1989-01-08 13:44:11.06 
17 65.771 20.312 0.09315 0.01804 1989-01-08 22:51:55.90 
18 65.778 20.314 0.09317 0.01802 1989-01-08 22:57:00.89 
19 65.824 20.322 0.09343 0.01783 1989-01-08 23:26:56.28 
20 67.294 20.603 0.08905 0.01612 1989-01-09 15:12:26.52 
21 67.759 20.682 0.08797 0.01660 1989-01-09 20:12:27.10 
22 67.843 20.698 0.08763 0.01633 1989-01-09 21:09:52.73 
1 49.879 16.466 0.13272 0.03920 1989-01-03 05:12:11.54 
2 49.895 16.470 0.13272 0.03914 1989-01-03 05:19:11.45 
3 54.889 17.836 0.12134 0.03259 1989-01-04 19:30:56.18 
4 55.010 17.869 0.12098 0.03216 1989-01-04 20:31:05.98 
5 55.056 17.881 0.12090 0.03200 1989-01-04 20:53:56.28 
6 55.177 17.913 0.12085 0.03160 1989-01-04 21:53:56.57 
7 58.566 18.754 0.11190 0.02726 1989-01-06 02:14:55.70 
8 58.964 18.849 0.11252 0.02601 1989-01-06 05:48:08.09 
5 9 60.786 19.268 0.10635 0.02415 1989-01-06 22:04:55.80 
10 60.812 19.274 0.10636 0.02406 1989-01-06 22:19:56.09 
11 61.300 19.383 0.10476 0.02446 1989-01-07 02:42:56.18 
12 61.317 19.387 0.10464 0.02439 1989-01-07 02:52:55.80 
13 61.608 19.453 0.10353 0.02324 1989-01-07 05:40:25.85 
14 61.625 19.457 0.10352 0.02318 1989-01-07 05:50:26.33 
15 64.878 20.139 0.09491 0.01917 1989-01-08 13:37:25.85 
16 64.888 20.141 0.09492 0.01913 1989-01-08 13:44:11.06 
17 65.771 20.312 0.09315 0.01804 1989-01-08 22:51:55.90 
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表 3 续 
Table 3 Continued 
Arc Number a/® 6/° àa :st) 6/((") 71) Date (TDB) 
18 65.778 20.314 0.09317 0.01802 1989-01-08 22:57:00.89 
19 65.824 20.322 0.09343 0.01783 1989-01-08 23:26:56.28 
20 67.294 20.603 0.08905 0.01612 1989-01-09 15:12:26.52 
21 67.759 20.682 0.08797 0.01660 1989-01-09 20:12:27.10 
22 67.843 20.698 0.08763 0.01633 1989-01-09 21:09:52.73 
23 67.929 20.716 0.08717 0.01579 1989-01-09 22:10:56.09 
24 67.945 20.717 0.08741 0.01599 1989-01-09 22:19:56.95 
25 68.014 20.732 0.08731 0.01551 1989-01-09 23:08:56.28 
26 69.727 21.027 0.08309 0.01482 1989-01-10 18:41:55.99 
27 69.839 21.047 0.08230 0.01449 1989-01-10 20:03:15.00 
28 69.872 21.053 0.08214 0.01438 1989-01-10 20:27:36.02 
29 69.876 21.053 0.08221 0.01437 1989-01-10 20:30:26.23 
30 70.094 21.090 0.08175 0.01367 1989-01-10 23:09:55.90 
31 70.115 21.094 0.08164 0.01456 1989-01-10 23:17:30.36 
32 70.176 21.104 0.08182 0.01343 1989-01-11 00:10:02.23 
5 33 70.705 21.192 0.07971 0.01317 1989-01-11 06:31:56.38 
34 71.893 21.380 0.07667 0.01230 1989-01-11 21:16:26.71 
35 71.908 21.383 0.07662 0.01225 1989-01-11 21:28:21.24 
36 73.016 21.567 0.07286 0.01112 1989-01-12 11:49:40.54 
37 73.025 21.568 0.07285 0.01108 1989-01-12 11:56:22.30 
38 73.027 21.554 0.07326 0.01123 1989-01-12 11:57:26.23 
39 73.039 21.558 0.07293 0.01126 1989-01-12 12:06:43.51 
40 73.058 21.561 0.07288 0.01117 1989-01-12 12:22:26.14 
41 73.087 21.565 0.07282 0.01105 1989-01-12 12:45:56.18 
42 73.093 21.566 0.07281 0.01102 1989-01-12 12:50:55.99 
43 73.094 21.566 0.07287 0.01098 1989-01-12 12:51:56.47 
44 73.097 21.565 0.07315 0.01095 1989-01-12 12:54:56.18 
45 73.166 21.577 0.07290 0.01069 1989-01-12 13:51:26.52 
46 73.224 21.585 0.07305 0.01047 1989-01-12 14:38:55.99 
47 73.767 21.661 0.07172 0.01051 1989-01-12 21:52:08.57 
48 74.739 21.798 0.06885 0.00994 1989-01-13 11:14:26.33 
49 74.844 21.813 0.06860 0.00948 1989-01-13 12:45:56.18 
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X) RAB: 基于 智能 优化 算法 的 小 天 体 初 轨 确 定 
表 3 续 
Table 3 Continued 
Arc Number a/® 6/° àf) sE) ÈC +877) Date (TDB) 
50 74.856 21.815 0.06859 0.00943 1989-01-13 12:56:56.28 
51 79.426 22.382 0.05673 0.00637 1989-01-16 11:22:33.62 
52 79.427 22.383 0.05673 0.00636 1989-01-16 11:24:03.48 
53 79.428 22.383 0.05673 0.00635 1989-01-16 11:25:35.93 
54 79.563 22.387 0.05667 0.00592 1989-01-16 13:46:11.16 
55 79.571 22.387 0.05668 0.00589 1989-01-16 13:54:26.23 
i 56 79.578 22.388 0.05670 0.00586 1989-01-16 14:02:41.30 
57 79.623 22.393 0.05686 0.00568 1989-01-16 14:49:25.85 
58 79.625 22.393 0.05687 0.00567 1989-01-16 14:51:56.18 
59 80.801 22.522 0.05350 0.00554 1989-01-17 11:15:36.31 
60 80.805 22.523 0.05348 0.00552 1989-01-17 11:20:36.12 
61 80.810 22.523 0.05347 0.00549 1989-01-17 11:25:35.93 
62 87.926 23.036 0.03938 0.00243 1989-01-23 17:21:52.15 
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